闲话/社论 22.12.15 两道计数题

闲话

好今天闲话终于有点拿得出手的东西了
向计数巨佬们致以敬意和谢意。

既然代码之后补了，那我这也少写点闲话吧（
一会再补！

下面是社论部分！

两道计数题

首先感谢 jijidawang 和 Alpha1022。

以及还有人记得这题吗？22.10.05写的那道！

CF1515E

给定 $n$，$M$。你有 $n$ 台电脑排成一排，你需要依次开启所有电脑。

你可以手动开启一台电脑。在任意时刻，若电脑 $i-1$ 与电脑 $i+1$ 都已经开启 $(1<i<n)$，电脑 $i$ 将立刻被自动开启。你不能再开启已经开启的电脑。

求你有多少种开启电脑的方案。两个方案不同当且仅当你手动开启的电脑的集合不同，或是手动开启电脑的顺序不同。答案对 $M$ 取模。

我们首先枚举有 $k$ 个位置的电脑是自动开启的，并将每一段全手动开启的电脑的方案数用 egf 表示。
容易发现答案可以被表为

$$\sum_{k=0}^n \left[\frac{x^{n-k+1}}{(n - k + 1)!}\right] \left(\sum_{i\ge 1} \frac {2^{i-1}}{i!} x^i \right)^k = \sum_{k=0}^n (n - k + 1)! [x^{n - k + 1}] \left(\frac{e^{2x} - 1}{2} \right)^k$$

考虑设

$$f = \sqrt{\frac{x(e^{2x} - 1)}{2}}$$

则 $f$ 的最低次项为 $x$，系数为 $1$。对于一个 $k$，我们要计算的就是 $[x^{n+1}] f^{2k}$。考虑拉格朗日反演，设 $g = f^{\langle -1\rangle}$，则我们要求的就是

$$[x^{n+1}] f^{2k} = \frac {2k}{n + 1} [x^{n - 2k + 1}]\left(\frac{x}{g}\right)^{n+1}$$

如果我们能求出 $g$，我们就可以迅速得到每个值。

由 $g$ 的定义得到

$$\sqrt{\frac{g(e^{2g} - 1)}{2}} = x$$

即

$$\frac{g(e^{2g} - 1)}{2} - x^2 = 0$$

分治 fft / 多项式牛顿迭代即可。总时间复杂度 $O(n \log n)/O(\frac{n \log^2 n}{\log \log n})$。

Submission.

牛迭

poly newton(int n) {
    poly F(2), A, B; F[1] = 1;
    for (int k = 2; k < (n << 1); k <<= 1) {
		F.redegree(k);
		A = (2 * F).exp();
		B = ginv(2) * (A - 1) + F * A;
		A = ginv(2) * ((A - 1) * F) - "x^2"_p;
		while(B[0] == 0) B <<= 1, A <<= 1;
		F -= A * B.inv();
    } return F.split(n);
}

signed main() {
    cin >> n;
	poly g = newton(n + 2);
	g = (g << 1).inv().pow(n + 1);
	for (int k = 0; n - 2 * k + 1 >= 0; ++ k) {
		ans = (ans + 1ll * gfac(n - k + 1) * 2 * k % mod * ginv(n + 1) % mod * g[n - 2 * k + 1]) % mod;
	} cout << ans << '\n';
}

---

[PA 2019 Final] Grafy

计数 $n$ 个点的有标号无向简单图，满足每个点的入度、出度均为 $2$，且无重边、自环。

$n \le 500$。

不妨设标号为 $0\sim n - 1$。

考虑计数边集组成的排列。我们将每条边拆成两条无向边，随后抽象出边为长 $2n$ 的排列 $0\sim 2n - 1$。$p_u = v$ 表示 $\left\lfloor u/2 \right\rfloor \to \left\lfloor v/2 \right\rfloor$ 有边。
一张图同一点的两条出边可以换位置，表示 $v$ 节点的值（$2v, 2v + 1$）也可以换位置，这样每张图能导出 $2^{2n}$ 个排列。

怎么算呢？考虑一手容斥。

思考什么样的排列是不满足条件的。

1. 重边。
   $\left\lfloor {p_{2i}}/2 \right\rfloor = \left\lfloor p_{2i + 1}/2 \right\rfloor$ 肯定不满足。
2. 自环。
   $\left\lfloor {p_{2i}}/2 \right\rfloor = i$ 或者 $\left\lfloor p_{2i + 1}/2 \right\rfloor = i$ 都不可以。

因此答案容斥信息可以表示为（这不是很懂）：

$$\prod_{i=0}^{n-1} \left(1 - [重边] - [2i\ 自环] - [2i+1\ 自环] + 2[均自环]\right)$$

随后可以枚举不满足条件的点的数量，即出现了 $i$ 个双自环，$j$ 个单自环，$k$ 个重边。有答案

$$\small\frac{1}{2^{2n}} \sum_{i + j + k \le n} \binom{n}{i} \binom{n-i}{j} \binom{n-i-j}{k} 2^i(-1)^{j+k}\times 2^i\times 2^{2j} \times (n - i - j)^{\underline k} 2^k \times (2n - 2i - j - 2k)!$$

四个 $\times$ 分割的五个部分分别是：容斥系数；双自环的贡献；单自环的贡献；重边的贡献；剩下部分的贡献。

$$\frac {1}{2^{2n}}\sum_{i + j + k \le n} (-1)^{j+k} 2^{2i+2j+k} \frac{n!\times (n - i - j)! \times (n - i - j - k)!}{i!j!k!(n-i-j-k)!^2}$$

枚举计算即可。复杂度 $O(n^3)$。

code

全屏放得开（

int n, ans, pw2[N];

signed main() {
	cin >> n >> mod; Mod.init(mod);
	init_fac(n << 1);
	pw2[0] = 1; rep(i,1,n<<1) pw2[i] = add(pw2[i - 1], pw2[i - 1]);
	rep(i,0,n) rep(j,0,n-i) rep(k,0,n-i-j) {
		addi(ans, mul(fac[(n - i - k) * 2 - j], ifac[i], ifac[j], ifac[k], ifac[n - i - j - k], pw2[(i + j) * 2 + k], (j + k) & 1 ? mod - 1 : 1, fac[n - i - j], ifac[n - i - j - k]) );
	} muli(ans, fac[n]);
	muli(ans, qp(pw2[n << 1], mod - 2)); // act as a prime
	cout << ans << '\n';
}

能不能再给力点啊？

$n \le 10^5$。

换种统计答案的方法。

我们由一张原图建一张新图，新图是一张有重边无自环的边带标号图。假设原图内有 $(i,u), (i,v)$ 两条边，那么在新图中连一条 $(u,v)$ 边，标号为 $i$。由于边无法退化成点，因此保证了原图无重边，这样只需要原图任意边的标号不是自己的一个端点，且没有相同标号就行了。再观察我们可以发现新图是一系列环拼合起来的。自然考察单个环。

我们仍然沿用容斥的思路，现在就需要钦定一系列边的标号是自己的一个端点。考虑令 $G(x, y)$ 为一个环的 egf，令 $[x^n y^k / n!] G(x, y)$ 为 $n$ 个点的环的计数，满足有 $k$ 条边未被确定标号。
先不考虑未确定标号的边，留下来的就是一系列链了。可以发现一条链上一定有一个点满足其左侧边需要选左端点，右侧边需要选右端点。因此假如有一条 $i$ 条边 / $i + 1$ 个点的链，其对答案的贡献是 $i + 1$。然后环可以翻转所以除 $2$，各段也可以顺次转动所以乘 $n / k$。
所以这部分（$k > 0$）的系数就是

$$\frac {n}{2k}{{n}\brack{1}} [x^n]\left(\sum_{i>0} ix^i\right)^k = \frac{n!}{2k} [x^n]\left(\frac{x}{(1 - x)^2}\right)^k$$

$k = 0$ 时平凡，就是第一类斯特林数。

于是

$$\begin{aligned} G(x, y) & \ = \sum_{n \ge 2}\sum_{k \ge 1} \frac{x^n y^k}{n!} \frac{n!}{2k} [x^n]\left(\frac{x}{(1 - x)^2}\right)^k + \sum_{n\ge 2}\frac{x^n}{n!} (n-1)! \\ & \ = \frac 12 \sum_{n \ge 2}\sum_{k \ge 1} \frac{x^n y^k}{k} [x^n]\left(\frac{x}{(1 - x)^2}\right)^k + \sum_{n\ge 2}\frac{x^n}{n} \\ & \ = \frac 12\left( \sum_{k \ge 1} \frac{y^k}{k}\left(\frac{x}{(1 - x)^2}\right)^k - xy\right) + \sum_{n\ge 2}\frac{x^n}{n} \\ & \ = \frac 12\left( \sum_{k \ge 1} \frac{\left(xy / (1 - x)^2\right)^k}{k} - xy\right) + \ln \frac 1 {1 - x} - x \\ & \ = \frac 12\left( \ln \frac 1{1 - xy / (1 - x)^2} - xy\right) + \ln \frac 1 {1 - x} - x \end{aligned}$$

这是一个环的 egf。得到原图的 egf $F(x,y)$ 只需要将它作 $\exp$，即

$$\begin{aligned} F(x, y) & \ = \exp G(x, y) \\ & \ = \exp\left(\frac 12\left( \ln \frac 1{1 - xy / (1 - x)^2} - xy\right) + \ln \frac 1 {1 - x} - x\right) \\ & \ = \frac{\exp(-\frac 12 xy - x)}{(1 - x)\sqrt{1 - xy / (1 - x)^2} } \end{aligned}$$

经典枚举标号随意的边数为 $i$，容斥系数为 $(-1)^{n - i}$。答案即为

$$\begin{aligned} & \sum_{k=0}^n (-1)^{n - k} \left[\frac{x^ny^k}{n!k!}\right] F(x, y) \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {(-1)^{n - k} n!k!}\right] F(x, y) \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] F(-x, -y) \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{\exp(-\frac 12 xy + x)}{(1 + x)\sqrt{1 - xy / (1 + x)^2} } \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \exp(-xy / 2) \times \left({1 - xy / (1 + x)^2}\right) ^ {-1/2} \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0}\frac {(-xy / 2)^i}{i!} \times \sum_{j\ge 0} \binom{-1/2}{j} \left(- \frac {xy}{(1 + x)^2}\right) ^ {j} \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0} \sum_{j\ge 0} \frac {(-xy / 2)^i}{i!} \frac{(-1/2)^{\underline j}}{j} \left(- \frac {xy}{(1 + x)^2}\right) ^ {j} \\ =\ & \sum_{k=0}^n \left[\frac{x^ny^k} {n!k!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0} \sum_{j\ge 0} y^{i+j}\frac {(-x / 2)^i}{i!}\frac{(-1/2)^{\underline j}}{j!} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{i\ge 0} \sum_{j\ge 0} (i + j)! \frac {(-x / 2)^i}{i!}\frac{(-1/2)^{\underline j}}{j!} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} & (k = i + j) \\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \sum_{i\ge 0} \frac{(i + j)!}{j!} \frac {(-x / 2)^i}{i!} & (提出\ i\ 相关部分) \\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \sum_{i\ge 0} \binom{i+j}{i}(-x / 2)^i \\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} \sum_{i\ge 0} (-1)^i \binom{-j-1}{i} (-x / 2)^i & (上指标反转) \\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2}\right) ^ {j} (1 + x / 2)^{-j-1} \\ =\ & \left[\frac{x^n}{n!}\right] \frac{e^x}{(1 + x)(1 + x / 2)} \times \sum_{j\ge 0} (-1/2)^{\underline j} \left(- \frac {x}{(1 + x)^2(1 + x / 2)}\right) ^ {j} \end{aligned}$$

气势磅礴……来自 EI。这里的推导尽量没跳步。

最后的式子里唯一不像微分有限的是那个 $(-1/2)^{\underline j}$。然而需要了解的是，

$$P(x) = \sum_{j\ge 0} (-1/2)^{\underline j} (-x)^j = \sum_{j\ge 0} (-1)^j (-1/2)^{\underline j} x^j = \sum_{j\ge 0} (1/2)^{\overline j} x^j = F\left( \left.\begin{aligned} 1, &1, \frac 12 \\ &1 \end{aligned} \right\rvert x \right)$$

超几何函数形式。其满足的递推形式为

$$x^2 P'(x) = (1 - x / 2) P(x) - 1$$

因此最后的整个式子都是微分有限的。根据经典结论，递推可做到 $O(n) / O(\sqrt n \log n)$。

code

#include <bits/stdc++.h>
using namespace std; using pii = pair<int,int>; using vi = vector<int>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
template <typename T> T rand(T l, T r) { return uniform_int_distribution<T>(l, r)(rnd); }
template <typename T> void get(T & x) {
	x = 0; char ch = getchar(); bool f = false; while (ch < '0' or ch > '9') f = f or ch == '-', ch = getchar();
	while ('0' <= ch and ch <= '9') x = (x << 1) + (x << 3) + ch - '0', ch = getchar(); f && (x = -x); 
} template <typename T, typename ... Args> void get(T & a, Args & ... b) { get(a); get(b...); }
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
template <typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
template <typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define ufile(x) 
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define Aster(i, s) for (int i = head[s], v; i; i = e[i].next)
#define all(s) s.begin(), s.end()
#define eb emplace_back
#define pb pop_back
#define em emplace
const int N = 1e7 + 5;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;

int mod, inv2;
// const int mod = 10007;
// const int mod = 469762049, g = 3;
// const int mod = 998244353; // const int g = 3;
// const int mod = 1004535809, g = 3;
// const int mod = 1e9 + 7;
// const int mod = 1e9 + 9;
// const int mod = 1e9 + 3579, bse = 131;
/* add / sub */ 			template <typename T1, typename T2> T1 add(T1 a, T2 b) { return (a += b) >= mod ? a - mod : a; } template <typename T1, typename ...Args> T1 add(T1 a, Args ... b) { return add(a, add(b...)); } template <typename T1, typename T2> T1 sub(T1 a, T2 b) { return (a -= b) < 0 ? a + mod : a; } template <typename T1, typename ...Args> T1 sub(T1 a, Args ... b) { return add(a, add(b...)); } template <typename T1, typename T2> void addi(T1 & a, T2 b) { (a += b) >= mod ? (a -= mod) : true; } template <typename T1, typename ...Args> void addi(T1 & a, Args ...b) { addi(a, add(b...)); } template <typename T1, typename T2> void subi(T1 & a, T2 b) { (a -= b) < 0 ? (a += mod) : true; } template <typename T1, typename ...Args> void subi(T1 & a, Args ...b) { subi(a, add(b...)); }
/* Fastmod / mul */ 		struct FastMod { int m; ll b; void init(int _m) { m = _m; if (m == 0) m = 1; b = ((lll)1<<64) / m; } FastMod(int _m) { init(_m); } int operator() (ll a) {ll q = ((lll)a * b) >> 64; a -= q * m; if (a >= m) a -= m; return a; } } Mod(mod); int mul(int a, int b) { return Mod(1ll * a * b); } template <typename T1, typename T2> int mul(T1 a, T2 b) { return Mod((long long)(1ll * a * b)); } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); }  template <typename T1, typename T2> void muli(T1 & a, T2 b) { a = Mod(1ll * a * b); } template <typename T1, typename ...Args> void muli(T1 & a, Args ...b) { muli(a, mul(b ...)); } // /* trivial multiple function(idk whether its useful*/ int mul(int a, int b) { return 1ll * a * b % mod; } template <typename T1, typename T2> int mul(T1 a, T2 b) { return (long long)(1ll * a * b) % mod; } template <typename T, typename ...Args> int mul(T a, Args ...b) { return mul(a, mul(b...)); }

int n, m, fac[N], inv[N], f[N], g[N], t1, ret;
int qp(int a, int b) {
	int ret = 1;
	while (b) {
		if (b & 1) ret = mul(ret, a);
		a = mul(a, a);
		b >>= 1;
	} return ret;
}

vi operator*(vi a, vi b) {
	static long long tem[100001];
	int size = a.size() + b.size() - 1;
	vi c(size);
	memset(tem, 0, size << 3);
	for (int i = 0; i < a.size(); i++) for (int j = 0; j < b.size(); j++)
		tem[i + j] = add(tem[i + j], mul(a[i], b[j]));
	for (int i = 0; i < size; i++) c[i] = tem[i] % mod;
	return c;
}

vi operator-(vi a, vi b) {
	if (a.size() < b.size()) a.resize(b.size());
	for (int i = 0; i < b.size(); i++) a[i] = sub(a[i], b[i]);
	return a;
}

vi deri(vi a) {
	for (int i = 0; i < a.size() - 1; i++) a[i] = mul(a[i + 1], i + 1);
	a.pop_back();
	return a;
}

vi shift(vi a, int n) {
	a.resize(a.size() + n);
	for (int i = a.size() - 1; i >= n; i--) a[i] = a[i - n];
	rep(i,0,n-1) a[i] = 0;
	return a;
}

void div(int a[], int v) {
	rep(i,1,n) a[i] = sub(a[i], mul(v, a[i - 1]));
}

int C(int n, int m) {
	if (n < m or m < 0) return 0;
	return mul(fac[n], inv[m], inv[n - m]);
}

int main() {
	get(n, mod); Mod.init(mod); inv2 = mod + 1 >> 1;
	fac[0] = fac[1] = inv[0] = inv[1] = 1;
	rep(i,2,n) fac[i] = mul(fac[i - 1], i), inv[i] = mul(mod - mod / i, inv[mod % i]);
	rep(i,2,n) inv[i] = mul(inv[i - 1], inv[i]);
	rep(i,0,n) g[i] = inv[i];
	div(g, 1), div(g, inv2);
	vi H = {1, 2 + inv2, 2, inv2}, a = H, b = vi() - H, c = (H - shift(deri(H), 1));
	a = shift(a, 2), b[1] = add(b[1], inv2);
	b = b * c, c = c * H;
	f[0] = 1;
	rep(i,1,n) {
		t1 = 0;
		for (int j = 0; j < a.size(); j++) if (i + 1 - j >= 0)
			t1 = add(t1, mul(i + 1 - j, f[i + 1 - j], a[j]));
		for (int j = 0; j < b.size(); j++) if (i - j >= 0)
			t1 = add(t1, mul(f[i - j], b[j]));
		if (c.size() > i) t1 = add(t1, c[i]);
		f[i] = mul(mod - qp(b[0], mod - 2), t1);
	}
	rep(i,0,n) ret = add(ret, mul(f[i], g[n - i]));
	printf("%d\n", mul(ret, fac[n]));
	// printf("%.6lf\n", 1. * clock() / CLOCKS_PER_SEC);
}